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We formulate the theory of electron transport through coupled-quantum dots by extending the 
auxiliary operator representation. By using the generating functional technique, we derive the exact 
expressions for currents, dot-occupation numbers and spin correlations, and examine them based on 
the non-equilibrium Green's function method under the non-crossing approximation (NCA). Our 
formulation generalizes the previous NCA approaches by allowing full occupation numbers with a 
finite Coulomb repulsion. 
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I. INTRODUCTION 

Transport properties of the double quantum-dot system have been extensively studied both experimentally and 
theoretically. [1| This artificial molecule, analogous to the two-impurity Anderson problem, provides a good platform for 
examining the exciting physics of the correlated electron behavior, such as the Kondo effect. The double quantum- 
dot structure is also a fascinating subject from the point of view of possible applications in quantum computation, 
where it is suggested as a basic building block, with qubits being represented by electron spins in each quantum 
dot.@. 

The rich variety of correlated electron phenomena of a double quantum-dot system emerges from the inclusion 
of the electron-electron interactions. frjl-fl2j In this regard, many theoretical approaches have been developed that 
concentrate on a limited range of relevant parameter values, such as infinite Coulomb repulsion, a finite interaction 
but under equilibrium transport conditions, or symmetric dot occupations. However, in order to describe the control 
and measurement of quantum bits in detail. [l3l [l4| it is necessary to develop the theory that deals with correlated 
electron behavior in a wide range of interaction parameters, level occupancies in dots, as well as to include the 
time-dependent perturbations. The need for such a theory comes from the fact that in experimental studies of spin 
blockade in lateral coupled quantum dots, independent tunnel barrier tuning with arbitrary dot occupations have been 
achieved. [TH [l6j Similarly, initialization and manipulation of quantum bits requires description of sudden changes of 
energy levels in the dots due to the time varying bias and gate voltages. 

In this paper, we derive the expressions for the current, densities of states, dot occupancies, and spin correlations of 
the double-dot system. Our approach enables the treatment of arbitrary Coulomb interactions, occupation numbers, 
finite temperature, as well as time varying voltages. To do this, we extend the auxiliary operator representation 
and apply the non-equilibrium Green's functions method associated with the generating functional described by the 
non-crossing approximation. [l7i [l8| In order to establish the validity of our approach, we compare our results, in a 
variety of situations, with the previous NCA and exact methods like the numerical renormalization group (NRG) 
scheme. We find that our formulation reproduces the previous NCA results, but deviates from the NRG method. 
This is not surprising since it is well known that the NCA fails in describing the low-energy Fermi-liquid regime. Since 
the vertex corrections cure the low temperature transport properties. fl9| this work may be used for more involved 
further studies. 

The paper is organized as follows: In Sec. II we introduce the Hamiltonian of the double-dot system, and reformulate 
it in terms of the auxiliary particle operators in order to calculate the non-equilibrium Green's function associated 
with the generating functional. We derive the expressions for physical quantities via the relevant projection in the 
auxiliary particle occupation number subspace. In Sec. Ill, the transport properties of the double-dot system are 
examined by using the numerical calculations, and are compared with the previous results. We summarize our main 
results in Sec. IV. Some mathematical details are deferred to the the appendices. 
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II. CALCULATION METHOD 
A. Hamiltonian 

We model the system, consisting of two quantum dots connected in series to the left and right electrodes, by the 
Hamiltonian, 

n = n d ots + nieads + n T . (i) 

Taking the full electron-electron interaction into account, [13] the Hamiltonian of the coupled quantum dots is given 
by, 

T-Ldots = X/ | €a ° rLa ° + t HC i aa c & a + ytJ aff n a f + - ill i - — J ^ n a<7 n &a ' 1 - JS L ■ Sr, (2) 

aa \ ~ ^ ' a' ) 

where we assume that each dot (a = L,R) has energy levels e aa labeled with spin index (er =t, i), and is coherently 
coupled to the other (a) with the tunneling matrix element tn- The dot number operators are given in terms 
of the creation (annihilation) operators cJj CT (c Q(T ) by n ar7 = c^ aa c arT , and the spin operators are given by S a — 
\ So-cr' ^aa' c \ta c aa' : where a are the Pauli matrices. U a and Ui are the Coulomb interaction parameters for electrons 
on the dot a and inter-dot, respectively, while J is the exchange coupling constant. 

The second and the third term in Eq. (p} describe the Hamiltonian of the leads and tunneling between the dots 
and the adjacent leads; 

Pleads ^ ^ tkaa&kafj&kaa 7 

Kt = { T ka a l a v C <x<* + r fc? c Lofc£»ff} » ( 3 ) 

where ot„_ (o-kaa) creates (annihilates) an electron at the lead a, and the constants T£ a provide the coupling strength 
between the dot and the adjacent lead. We assume that the energy levels of each dot are controlled independently 
by the nearby gate electrodes[l| and the chemical potential /i Q at the lead a is adjusted by applied voltage difference 
A/it to be 

Mi = Mo - A/i/2, 

fi R = mo + A/V2, (4) 

with respect to the equilibrium chemical potential (1q. 

In order to take into account all the possible occupancies of the double-dot system (ranging from zero to four) , we 
extend the idea of the auxiliary particle representation (l9l l2ll] and introduce the auxiliary operators dj„(d m ) as, 

C L = Cmm'dlndm' ■ (5) 
mm' 

Here, |vac) is chosen as the m— th basis vector diagonalizing the isolated coupled-quantum dots (or molecular 
states) as specified in Table HI The auxiliary operators satisfy the commutation relations [d m ',dj n ] s = 5 m ,m' and 
[d m >,d m ]<; — [d m i, djrjs = 0, where odd (even)-numbered states are assumed to be fermionic (bosonic) and ? = + if 
both m and m! denote fermionic states, otherwise ? = — . The overlap matrix t, m a m t = (m \ c^ a(J \ m!) has a non-zero 
value for the combination of boson-fermion or fermion-boson operators to ensure original commutation relation of 
[c-aai c a 'a'] = Ga,a' 5a,<r' under the constraint of Q = J^m^m^m = 1- The proof is given in Appendix lAl In terms of 
the auxiliary operators the Hamiltonian of the coupled dots is given by, 

15 

Hdots = emd m d m + Hint j{d m ', din}). (6) 

m=0 

The term Hint represents the interaction between auxiliary particles (applicable for Q > 2) and we omit it hereafter 
since it doesn't affect our final results. 
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Table I: Schematic representation of the eigenstates for two coupled quantum dots with a coupling strength tn- Here, energy 
splitting for one- and three-particle states are given by Ai CT = -J (ejj CT — tLa) 2 /4+ify and A3o- = \J (eLa + UL — £n(T — Un) 2 /^+t 2 H , 
respectively. Energies of two-particle singlet states are solutions of a cubic equation, Ae 3 +(UL+Un,— Ui — J)(Ae 2 — [El — En] 2 ) — 
([El — Er] 2 + 16t 2 H )Ae = and their eigenstates are determined by u k = (Ae k + El — En)w k /y/St H and v k = (Aek + UL+Un — 



Ui — J)uk/V&tH — Wk under a normalization of uf. + ■ 



— 1. We also abbreviate E a = e a f + e a ^ + U a 



0- particle state 

1- particle state 



2-particle state 



3- particle state 

4- particle state 



jm = 0) =| e) 

I "1 = 1,2} = (cos cj> a c La -sin cj> a c Rcl ) 
|m = 3,4) = (sin^c^+cos^c^) 
with tan20 CT = 2t H /{^n<r — £l<t) 
|m = 5> = cl t c R1 . | e) 
\m = Q) =c^4i I e ) 

i m = 7 ) = 72( c it c fli+ c L4 t ) i e ) 

where 5 S = ^(cj^c^-c^c^), 
jm = ll,12} = (cos9 a c L <T— smOaCna) \ f) 
jm = 13, 14} = (swOvClv+cosOvCrv) \ f) 
with tan20 CT = 2t H / {^na + Un — £l<t — Ul) 
| m = 15} = 4 t 4 4 c}j t 4 i | e} =|/} 



eo = 

£l,2 = (e£, CT +efl CT )/2-Ai CT 
£3,4 = (eL CT + efl CT )/2 + Ai CT 

e 5 = £ L t + £flT + (^- J )/ 2 
£ 6 = £L4. + £fl| + ([//-J)/2 
£7 = (e 5 + £6)/2 
£fc = ( J B L + £ fl )/2 + Ae fe /2, 
fc = 8,9, 10 

£11,12 = (£15+£Lct+£Hct)/2— A3 CT 
£13,14 = (£l5+£L t T+£H<r)/2 + A8 CT 

£15 = £L+£fl + 2(//-J 



B. Generating functional 

For the ease in evaluating the expectation value of any operator O, we introduce a Lagrange multiplier A associated 
with the auxiliary particle number Q as, 

n^n + xQ-Y, ^^ ads , (7) 

a 

with n l ^ ads = J2ka a kaa a kaa- Then, the system becomes the grand canonical ensemble with respect to the auxiliary 
particle number Q, i.e., Q is now unconstrained. 

With the grand canonical ensemble, we define a generating functional W = — InZ as an extension of the Gibbs 
free energy. Here, the generalized partition function Z, in terms of the coherent path integral representation is given 

Z = j> V[c* aa , c a<y ,a* kaa ,a kaa ]e~ s / lh , (8) 
with the action represented on a closed time contour as, 

S = j> dr [" d* m (t) (ihdT -e m -\) d m (t) + ^ a* ka(J (r) (iWr - e toff + /j, a ) a ka(7 (r) 

- E { T *>L>)<w(r) + T£c* aa (T)a ka<r (r)}] . (9) 

We note that the Fermi (Bose) particle operators are now replaced by the corresponding Grassman (complex) variables 

{ a acri a aai d m , dm}. 

From the unconstrained generating functional W, the expectation value in the Q — 1 ensemble can be calculated 
easily by noting that the operator Q commutes with the total Hamiltonian, and Q is thus a good quantum number. 
This fact enables us to expand the partition function in powers of £ = e~ A/9 . The partition function belonging to 
the Q = 1 subspace can be obtained by differentiating with respect to £, i.e., Zq = \ = limx->-oo g^-Z- Based on this 
relation, we can evaluate the expectation value of O by taking a functional derivative of Z with respect to its conjugate 
variable rj as, 



(0)c = -— !— ^-Z Q =x = lim 
£q=i or) A->-oo 



, (d/dO(o) GC 

\y)Gc - 



(10) 



where (0)c denotes the average over the canonical ensemble, i.e., over the subspace Q = 1 while (O)gc = SW/Srj is 
the average over the grand canonical ensemble. When the canonical expectation value of the operator O has a zero 
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expectation value in the Q = subspace, the above relation is further simplified to, 

(0) GC 



(0) c = Inn 



A-J-oo (Q)gC 



(11) 



As seen in the following section, since the expectation values of interest have a zero expectation value in the Q = 
subspace we hereafter focus on the average over the grand canonical ensemble based on Eq. (fTTj) . 

The partition function in the grand canonical ensemble is calculated following the standard series expansion 
procedure. [23j Firstly, we integrate the action over Grassman variables {a>% aa ,a-kaa} and obtain 



S = I drJ2d* m (T)(ihdT - e m -\)d m ( T )- idrdr 1 \ | 2 (t) 5jW (r, r')c aa (r'). 

Then we express the variables {c* aa ,c aa \ in Eq. ([5]) in terms of their functional derivatives as, 

8 8 



(12) 



in 111 ' 



acr* 

mm' ^m' 



Sv,u(t) Sri* m ,(T) 

8 8 
8r)m'(r) 8t]^(t) 



(13) 



where rj m and r\* m are their corresponding conjugate variables. In that way the partition function is further rewritten 

as, 



Z-Z°exp|l^ ^drdr'^^^^^rOwCToje^E™/^'^^^.')^^ 



(14) 



where Zq is the unperturbed partition function and <Jw(t, r') = | T? a | 2 gkaa( T iT~'). Here, the lower case g's are 
the unperturbed Green's functions of the lead-electrons and auxiliary particles; for instance, retarded, Keldysh, and 
advanced components of the Green's functions of the m-th auxiliary particle are given by 



(*,*') = ^(t-t'){[d m (t),dl(t')U m )° GC = Ifl( t 



ifl 

1 



= g**(t',t) 



- X[ d m (t),dl(t')},J GC = - 
in in 



tanh 



iti 

P(e m + A) 



a (e m +A)(i-t')/ifi 



(15) 



respectively, where the superscript '0' denotes the average in the case of the dots being decoupled from the leads, 
and <^ m is — 1(+1) if the particle to is a fermion (boson). To simplify the expressions, it is sometimes convenient to 
use greater, g > = (g K +g R — g )/2, lesser, g < = (g K — g R +g A )/2, and correlated, g c = g R — g A , Green's functions 
interchangeably 

Next, we expand the exponential function in Eq. (|14[) in power series and obtain the partition function by collecting 
all the connected diagrams; 



2 _ 2° e -{ w<1) + vv(2) + w(3> +'--} 
Here, is the collection of the | T£ |™-order diagrams. For instance, W^ 1 ' and W*- 2 -* look like, 



(16) 



ih 



W ( 2 ) = 




(17) 



where the solid (dotted) lines denote the unperturbed Green's functions g rn {g aa ). Here, large dots indicate the times 
at which the tunneling events occur and the overlap matrix is assumed to be multiplied as, 



; oca 



: ceo 



xC 



m'mi 



mm f 



(18) 
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Finally, the calculation of the generating functional W = — In Z is done in a straightforward fashion via the Luttinger- 
Ward functional $ and repeated terms; [24| 

W = -hxZ(°) + W (1) + W (2) + • • • 

= $ + V] <r p Tr [in Gp 1 + S P G P ] , p = acr and to. (19) 
p 

Here, the Luttinger-Ward functional $ is the sum of all the closed skeleton diagrams with a non-interacting Green's 
functions (g p ) replaced by the full Green's functions (G p ). 

Up to now, the generating functional of the coupled-dot system has been derived without any approximations, and 
thus the associated Green's functions give the exact expressions for the physical quantities as shown in the Appendices 
IB1 and [Cl In the next section we describe the approximation of the Luttinger-Ward functional <J>, and present the 
expressions for the physical quantities in the static case. 

C. Non-crossing approximation and projection to Q = 1 

Hereafter we employ the non-crossing approximation, that is, we confine our attention to the first skeleton diagram 
originated from W^ 1 ', and approximate the Luttinger-Ward functional by, 

m 

= -iTl ^™l&<| 2 j> ' dTdT'G m {T\T)G aa (T,T')G m ,(T,T') (20) 

mm'aer 

where the thick lines represent the full Green's functions of the particles instead of unperturbed ones (thin lines) in 
Eq. (H7D- 

Then, since the generating function W is stationary with respect to G p , namely 8W/8G P = 0, the self-energies can 
be obtained from, 

^ T ^ = -^5G^ry (21) 

Using the NCA functional Eq. (|2"0")) , this gives, 

S Q ,(r,r') = -ih^Sm |CW I 2 G m (T,r')G m ,(r',T), 

mm' 

E m (r,r') = ih [ ICV I 2 G QCT (r,r')- \^ >m \ 2 G aa (r> \r)]G m ,(r,T'), (22) 

m' aa 

for the self-energies of electrons in the leads and auxiliary particles. As functions of real time arguments, the above 
expressions can be rewritten as, 

£*(t,i';A) = -^E^l&'l 2 [G^(t,t';X)G<,(t',t-,X) + G<(t,t'-,X)Gi,(t',t;X)], 

mm' 

£*(M'; A) = -ih J2 «m ICV I' [ G ™(t,t'; A)G£,(t',i; A) + G°(t,t; X)G<,{t',t; A)] (23) 

mm' 

for the electrons in the leads, and 

£*(M';A) = ihJ2 [ICVfC(^;^)- \C^ m \ 2 G<At',t-X)]Gf nl (t,t';X), 

m ' a a 

S<(t,i';A) = ih [\Cn a m >\ 2 G<,(tJ;\)- |££J 2 ^(f.f; A)] G<,(t,t';X). (24) 

m'a<7 

for the auxiliary particles. 
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Next we project the self-energy to the Q — 1 ensemble (details can be found elsewhere (171 [la . |25|). For this we 
exploit two facts about the lesser and greater components of the Green's functions for the auxiliary particles. Firstly, 
it is important to note that self-energies of electrons in the leads depend on lesser components of Green function C?< . 
This makes the projection to Q = 1 subspace easy because 

G£(t,f) = ^(d m (t')U m (t)) GC = 0(e-^) (25) 
in 

means a zero expectation value in the Q = subspace, and thus one can use Eq. (|lip in evaluating observables. 
Secondly, the lesser and greater components of the Green's functions for the auxiliary particles are given by, 

G<'>(f,f) = /dtidt 2 G«(t,t 1 )E<' > (t 1 ,t 2 )G^(i 2 ,t') (26) 



without dependence on g^ > due to the loss of memory. (171 [20j| 

Before the projection to the Q = 1 subspace, we eliminate the A-dependence in g^ A . Since A is related to g^ A {t 1 1') 
only through the factor e A ^*~* in Eq. (fTijj) . the elimination of A such as g m (t,t';\) —> g m (t,t') = g m (t, t'; A = 0) 
results in the modified forms of self-energies in Eqs. $Z5§ and (JUJ) like, E p (t,t';X) t'\ \)e~ x( - t - t ">/ lh . 

By using Eq. (fTTj) and taking the projection of A — >• oo, the projected self-energies of Eq. (l23t become, 

E* (t.iO = ^E^l&'l 2 [G^(t,OG<,(t',i) + G<(t,OG^(t',0], 

mm' 

E£(t,f) = ^^ ?ra |^| 2 [G<(M')C(i'<) + G^(M')G<-(i',i)] (27) 

mm' 

where we use the abbreviated notation of 

V%?W)= lim e^S^^t's^e-^*-*')/^ 

A— voo 

and set a e^ x (Q)cc term aside in Eq. (1TT1) for a while. Here, G< (£,£') is defined in Eq. (|26| with its self-energy 
E<(i,t') ee lirru^oo e xfl T,<(t,t'; \) e -Mt-t')/^ _ Using Eq _ t h e se lf-energy is given by, 

E<(t,i') | 2 S< CT (MV IC',™I 2 G<,(t,0- (28) 

m'acr 

Whereas, the Dyson equation of G^ A (t,t') is 

G£ A (M') = <?£' A (M') + / dhdtzg^h^ih^G^itz,*) 
with its self-energy defined by E^ A (t,t') ee limA-**, E% A (t,t'; \) e -Ht-t')/th. from Eq ^ 

E£(t,i') = 0C%' I 2 <?^(M')- ICV I 2 <&(f,t)] G« (*,*')• (29) 

m'acr 

During the projection, we employ the relation, 

G< ff >(M';A) = g<>>{t,t') + / dt 1 dt 2 Gt(t,t 1 )E<^>(t 1 ,t 2] X)Gt(t2,t') (30) 



and neglect the second term, due to its 0(e dependence as seen from Eq. (|24|) . 
On the other hand, the expectation value of the operator Q is given by 

lim e A/5 (Q) GC = lim e A ^V <; m G<(t,t;X) = ihY* <; m G<{t,t) (31) 

A— >oo A— ¥oo z — ' ' — ' 

m m 

where the second step can be derived in a similar way to that of Appendix [Bj Throughout this work, we keep 
liniA^oo g x/3 (Q)gc to be unity via the normalization and consequently Eqs. (|27p - (j28|) are also the averaged values in 
the canonical ensemble. 

Eqs. (|27 p -(|28 p are the main results of this work, which can be applied to a double-dot system at arbitrary temper- 
ature, Coulomb interaction, source-drain and gate voltage configurations, including the time-dependent problems. 
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D. Physical quantities in static cases 

Since in the static case, the Green's functions depend only on the time interval, it becomes convenient to use the 
Fourier transform, 

1 r°° 

G(t,t') = — dEe E(t ~ t)/lh G(E). (32) 

By using the cut-off, p^ a (E) — 27r^ fc | Tg a | 2 S(E — eka), the unperturbed Green's function of electrons in the lead 
a, is then given in the energy space as, 

g< a >{E) = ±ihp^(E)f(±(E-^)), (33) 

where f(E) — 1/(1 + e@ E ) is the Fermi-Dirac distribution function. By defining the spectral function A m such that 
G^{t,t') — —2m<, m A m (t,t'), the self-energies of Eqs. (|27[l - (|2"8]l are rewritten as, 

/oo 
dE 1 [| Cn', m I 2 Gl,{E+E')- | | 2 GiAE'-E)} A m (E>), 

mm- -°° 

/OO 
dE' [| C',m I 2 Gg,(E+E')- | Z% m , | 2 Gg,(E'-E)] A m (E>) 

mm- -°° 

/oo 
dE' [i i 2 qIae-e')- i a, m i 2 g^-m] on m 

Here, the spectral function A Tn (E) is determined from 

/oo 
dE' [ICVI 2 g«*(E'-E)- I 2 

with the normalization condition of f dEA m (E) = 1 from Eq. (|3T|) . 

On the other hand, the expectation values of the physical quantities in the Q — 1 ensemble can be obtained by 
combining the results of the Appendix IBI with Eq. ((34]) . We summarize the results, in the energy space; for the 
current in the lead a, 

/oo 
dE [g> a (E)X< a (E) - g< a {E)T>AE)] , (35) 

for the density of states, 



for the occupation number, 



DOS aa {E) = --lmg^(E) = --lmZ% a (E), (36) 

7T 7T 



W)c = E ( IT*- ) / dEA m (E) , (37) 



dCouT J .1 x 



and for the spin-spin correlations, 



S, = & . s R)a = - s g£.^jwza r ^,„ (E) . (38) 



III. RESULTS AND DISCUSSION 



In this section, we illustrate the numerical solutions of Eq. (|34|) . and the resulting physical quantities of Eqs. ([35]) 
(|38|) , as well as the accuracy of the present theory. 
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For the cut-off function we choose a Lorentzian model p^ a (E) 

with W being the half width of the conduction band. 

In solving the Dyson's equations with self-energies given in Eq. (|34j) . we use the adaptive mesh scheme where more 
mesh points are inserted into a high weighted region for every interaction. The iteration is repeated until the following 
sum rules converge within 0.01%; 



1 f°° 

- / ImG*(E)dE = 1, 

7T J-oo 

/OO 
A m (E)dE = 1. 
-oo 



To achieve this numerical accuracy, we use about 1000 mesh points for each Green's function of an auxiliary particle. 

For simplicity we consider the symmetric case, Ul = Ur = U , tha = (r<j = (-d, Uj = 0, with T aa = T, and J = 0, 
and all the energies are measured in units of T (in an experiment V is typically of the order of fxeV to meV). We 
present results for two kinds of systems; one is a single quantum dot (that is, we take tji — > co) and the other is a 
double quantum dot (finite t#). Although the single quantum dot case has been extensively studied, we revisit the 
problem to show that our formulation indeed encompasses the previous results. 



A. Single quantum dot 

We first consider a single quantum dot and examine the correlated quantum transport through it. To do this 
we write td — >• £d + tn , and take tu — ► oo. Then, there are four low- lying states relevant to transport: | 0) =| e), 

I 1 ' 2 ) = 72^La ~ C Ra) I e )' and I 8 ) = Wl\ ~ C Rf)( C Ll ~ ^Ri) I e ) > Wnile their ener g ies are S iven b Y e = 0, £i, 2 = € d , 

and e 8 = 2e,j + \ {U L + U R + Uj + J), respectively. 

In the limit U —> oo, the state | 8) can be discarded further, and the present formalism recovers the results of Ref. 



17j. For this case, results of the system described by typical parameters are shown in Fig. [TJ In Fig. [T]-(a), we plot 
the equilibrium densities of states for several temperatures, where the broad peaks are caused by the usual transitions 
between levels (in this case | 0) and | 1, 2)), while the sharp ones (located at E = ^ = 2.0) are the Kondo peaks. The 
later ones increase as temperature is lowered. [l7j with saturation well below the Kondo temperature [lj| 

where 

1 1 



/ = HT 



| ea - li | ed + U - fi\ 



From this relation, Tk — 2.8 x 10~ 3 is estimated in the case of Fig. Q]while the calculated Kondo temperature, equal 
to its half width at half maximum, is 3.2 x 10~ 2 . 

The over-estimation of the Kondo temperature is a known consequence the NCA, as well as the Kondo peak 
height. [13] In the figure Fig. [T]-(b), we show the variation of the dot occupation, tie — (n m =i) + {n m =2) as a function 
of the dot energy level, and in Fig. E-(c), the relation between the electronic occupation ub and the height of the 
density of state at E = fi. Actually, ub and the density of state are related through the Friedel sum rule: 

DOS{E = no) = ^sin 2 (7m B /2). (41) 

7Tl 

For the validity of our calculations, we also plot previous results of the NCA and numerical renormalization group 
(NRG) method from Ref. [l8[ . In the case of the occupation, we find that our results are in good agreement with the 
previous results, with a minor deviation resulting from the use of a different cut-off function. On the other hand, in 
the comparison of the Friedel sum rule, a large deviation of our results are found from those of the exact result and 
NRG. As seen in results of Ref. [l8j], the previous NCA calculation also show nearly the same deviation. This fact 
leads us to the over-estimated Kondo peak with the NCA in the wider range of the occupation. 
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Figure 1: In (a), the equilibrium densities of states are plotted for a single quantum dot with an infinite repulsive potential. The 
inset shows the blow-up around Kondo peaks for temperatures T = 1.0 x 10 _4 (solid), 5.0 x 10 _4 (dotted), 1.0 x 10~ 3 (dashed), 
and 5.0 x 10 -3 (dot-dashed), respectively. We use the parameters of €d = —0.5, /io = 2.0, Afi — 0, and W = 100.0 which gives 
Tk = 2.8 x 10~ 3 in Eq. (|40p . In (b), we compare the electronic occupation ns as a function of energy level with the previous 
results. Crosses represent results from NRG, boxes from the NCA of Ref. [181 ] at T = and circles from our approach at 
T — 1.0 x 10 -4 , respectively. Since our calculation code is not available at T = 0, we choose sufficiently low temperature for 
the comparison. In (c), the sum rule is examined together with the exact result (solid line) from Eq. (I41[) . 



When the Coulomb interaction is finite, all of the four states | 0), 1 1), | 2), and 1 8) take part in electron transport. 
For the same quantum dot of Fig. [TJ but with a finite potential U, we show the DOS as a function of E in Fig. (J2J)-(a) . 
As the Coulomb potential decreases, the Kondo peak is lowered because the electron correlation is unimportant. This 
is also predicted by the Friedel sum rule of Eq. (|4ip . Since the two-particle state 1 8) becomes energetically favorable 
with smaller potential U, the electron occupation increases up to two. In Fig. H}(b), we plot the height of the Kondo 
peak as a function of occupation ub ■ It is found that the present result (circles) exhibits the same decreasing behavior 
with larger occupation as the Friedel sum rule, however still shows the over-estimation of the Kondo peak. In Fig. 
[2j(c), the spin correlation is shown as a function of the Coulomb potential U and is compared with the occupation 
ns- Since the spin correlation originates from only a two-particle state, results in the figures are proportional to the 
occupancy of 1 8) auxiliary particle. Thus, one can see that as the occupation tib approaches two, the spin correlation 
becomes S 2 = (8| S L ■ S R |8) = -3/8. 

Fig. [3] shows the linear conductance (a), and the electronic occupation tib (b) as a function of the single-particle 
energy for a finite Coulomb potential U = 10.0. For this calculation, we apply a small voltage between the left and right 
leads of A/x = 0.01 and the conductance is calculated as the current at a lead divided by A/i. As the temperature 
is lowered, the conductance increases and approaches 2e 2 //i, which is in accordance with the experimental results 
reflecting the Kondo effect. Q On the other hand, the electronic occupation ub shows weak temperature dependence, 
as shown in Fig. [3]-(b). The conductance maximum are approximately at tib = 0.5 and 1.5, which coincides with the 
condition of most probable sequential tunneling: i± = and \i + ex t 2 = eg. 



B. Coupled-quantum dots 



When the coupling strength tjj between the dots is finite, the system now represents double quantum dots and all 
16 many-body molecular states take part in the transport. As the first example, we consider the case of U — > oo and 
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Figure 2: In (a), we show the equilibrium density of states for the single quantum dot with the same parameters as Fig. [TJ 
but with finite potential U. In (b), the density of states at the chemical potential is plotted as a function of the occupation 
n_e, where the solid line is the Friedel sum rule of Eq. (|4ip and circles is the present result, respectively. In (c), we plot the 
variation of the occupation ns and spin correlation & with respect to the potential U. Dotted lines are the guide for eyes. 



Ui = J = 0. Then, low-lying states are eo = 0, ej.,2 = £d — tn, £3,4 = £d + iff, and e$fi,7,8 — 2e<i from Table HI Due 
to the large Coulomb potential U, one can see that the double occupation on each quantum dot is prohibited. And 
one expects that sequential tunneling occurs dominantly for two conditions of = tu and = —iff. The former 
corresponds to the transition between | 0) and | 1,2), and the latter is that between | 1, 2) and | 5, 6, 7, 8). 

In Fig. |4j we examine the conductance as a function of the chemical potential difference in the vicinity of the latter 
case. For a given = —2.5, we compare calculated conductance for — 2.0, 2.6, and 3.2. Among three cases, 
overall conductance for tn = 2.0 shows the largest value. It is interesting because the largest one will be the case 
tn = 2.6 according to the sequential tunneling condition of td = —iff- We attribute this to the level renormalization 
owing to the electron correlation. On the other hand, sharp peaks are found around Afj, = 0, whose height increases 
as temperatures are lowered. The peaks are found to result from the Kondo effect as inferred from the density of 
states in Fig. HJ-(b). 

Actually, the similar calculation is already performed in Ref. [l2j . where double peaks of the conductance around 
A/i = differently from the present result are observed. We attribute the discrepancy between a single peak and 
a double peak predicted in each work to the difference in the formulation of the problem. While in Ref. [12T |. the 
localized basis such as c^ aa | e) is used, we use the diagonalizing basis shown in Table fl] Strictly speaking, the present 
work treats the double-dot system as a single-quantum dot with multi-level molecular states, which leads to the 
modified coupling strengths between the dots and the leads, weighted by £mm' in Eq. (JS|). Therefore, although both 
approaches adopt the NCA, the details of the Feynman diagrams are different, and we expect that the results of both 
approaches would converge by including more crossing diagrams. 

Finally, we consider the coupled quantum dots with finite Coulomb potential. When the Coulomb potential becomes 
comparable to the dot-dot interaction tn , their competition gives rise to the rich electronic structure and all the energy 
levels may be relevant to the transport. In Fig. ©-(a), we show the variation of energy levels as a function of the 
Coulomb potential in the case of isolated coupled dots. As the Coulomb potential decreases from infinity, it is found 
that more levels fall into the range of relevant energy. In other words, this means that various transitions between 
states become available and are responsible for more peaks in the density of states as shown in Fig. ©-(b). Due 
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Figure 3: We plot the linear conductance (in unit of 2e 2 /h) as a function of single-particle energy in (a) and corresponding 
electronic occupation (b), for a finite Coulomb potential U = 10.0 at temperatures T = 0.003(solid), 0.03(dotted), and 



0.3(dashed), respectively. We use W 
difference A/i = 0.01 through I a /Afi. 



50.0 and (jlq = 4-0. The linear conductance is calculated with a finite potential 
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Figure 4: In (a), the conductance (in unit of 2e 2 /h) of the coupled-quantum dot with U = oo and Ui = are shown as a 
function of the chemical potential difference for temperatures T = 0.001 (thick), 0.005 (medium), and 0.01 (thin), respectively. 
In (b), we plot the density of states at A/x = 0.0 and T = 0.001. The value of the interdot interaction parameter for each of 
the curves is ta = 2.0 (solid), 2.6 (dotted), and 3.2 (dashed) in both panels and W = 10.0, = —2.5, and jito = 0. 



to the detailed change of energy levels, the conductance is also found to be largely modified. In Fig. ©-(c), the 
conductance are shown for three different Coulomb potentials. Compared to that of the infinite Coulomb potential 
case in Fig. (|l|)-(c), calculated results are largely suppressed. This is because the transition energies determined 
from the competition of various interactions are too large for electrons to tunnel through dots, which is similar to the 
Coulomb blockade effect for large potential U. Nevertheless, one can see sharp peaks in the calculated conductance 
at A/i = 0. These peaks result from the Kondo effect as in the case of the infinite Coulomb potential, meaning that 
correlated transport still occurs even in small Coulomb potential. We find that the heights of the conductance at 
Afi = are much larger than those of the master equation approach, however, smaller than those of NRG (not shown 
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Figure 5: For the coupled quantum dots of Fig. 0] the variation of energy levels is plotted in (a) as a function of Coulomb 
potential U where one-, two-, three-particle states are represented by solid, dotted, and dashed lines, respectively. The arrows 
indicate possible transitions into one-particle (£217*21') an< A three-particle (£23) state from the ground (two-particle single) state. 
In (b), we show the density of state for finite Coulomb potentials U = 3.0 (solid), 5.0 (dotted), and 10.0 (dashed), respectively, 
at temperature T = 0.003 and display the transitions corresponding to each peak. In (c), the conductance (in unit of 2e 2 /h) 
is shown as a function of the chemical potential difference. 

here) . [H [ll[ This means that our approach accounts for correlated behavior of electrons partially. 

IV. SUMMARY 



In summary, we formulate the electron transport through two laterally coupled quantum dots by extending the 
auxiliary operator method to a multi-level case, and derive the non-equilibrium Green's function in a conserving way. 
By using the generating functional technique, we present exact expressions for the current through the system, as well 
as the densities of states, occupancies, and spin correlations of the dots. To obtain the Luttinger-Ward functional, 
we include the first-order diagram (non-crossing approximation). For the validity of our results, we examine various 
situations and compare calculated results with those of previous NCA and exact NRG approaches. We find that our 
formulation encompasses the previous NCA results successfully, however gives the deviated behavior from the NRG 
method. This means that the present method accounts for the correlated behavior partially and the vertex correction 
is needed for more accurate description of transport. Nevertheless, since the present theory cope with all the ranges 
of the Coulomb energy and occupancies as well as time-dependent voltages, it can be applied to reveal transport 
properties of various double-dot problems. 
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Appendix A: commutation relations 



In this appendix we show that the auxiliary particle representation, 

C (Lcr = ^ ] £m,m'dmd m ' (Al) 
mm' 

gives the correct commutation relations [c aa ,c , ,]_)_ = S aa '6 aa i . Here, a quasi-particle operator d m is fermionic 
(bosonic) when the number of particles in a state it represents is odd (even), and is assumed to satisfy the commutation 
relation [d m ,d\ n ,]± = 5 mm i. It is important to note that the expansion coefficient S, m a m i = (m \ c\ a \ m') is nonzero 
only if the number of particles in | m) is larger than that in | m!) by one. This means that Eq. (|A1|) is the combination 
of fermion and boson operators. 

Now, we calculate the commutation relation, 

i c aa,C < a 'cr']+ = £,m 2 m 1 £m 3 m 4 [^m 2 ^ m i ' ^7713 ^'TiJ + • (A2) 

m\ ,7712 ,7713 >77l4 
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To calculate the right-handed side, it is convenient to separate the sums into the bosonic (B) and fermionic (F) terms, 

E = E (y, + Y,)(y, + Y,) 

">i,mj,ni3,m4 m 2 ,m 3 \mi£F roiSB/ \m 4 eF m 4 G-B/ 

e + eVe + e 

^mieF,m 2 eB m 1 eB,m 2 eF J \ m 4 e_F,m 3 e_B m4eB,m 3 eFi 

E E + E E V (A3) 

i mi,m3GF m 2 ,m^£iB m± ,m 3 £B 1112 ,m^GF J 

Here, in the second line we use the fact that c aa and cj^ are the products of fermionic and bosonic operators, and 
in the third line we display the collection of non-zero terms. By substituting Eq. (|A3j) into Eq. (|A2j) . and using 
[d m ,dl n ,]± = 6 mm >, we arrive at, 

\r = ( t a V f a '?! , 4- f ar7 ,* „f a ' a '„} d , (A4) 



mm m 



Furthermore, since the expansion coefficients satisfy the orthogonality relation, 

^ ] {^m"m^m"m' + £m' m" £mm") = ( m I [ c a<y 1 C q',t']+ I m ) — $mm' $ua' $aa' , (A5) 
m" 

the commutation relation is simplified to 

m 

Thus, in the subspace Q = 1, the combination of quasi-particle operators leads to the correct commutation relation 
between c aa and c^,^. 



Appendix B: Expressions for physical quantities 

In this section we derive the exact expressions for observables in terms of the Green's functions. As examples, we 
show the procedure for evaluating the current, 



M*) = ij t ( E 4« CT (*H-(*) ) = f E ( T £4r(*)<w(*) - T?A a Jt)c a At)) Gc , (Bl) 

\ k I GC k 



Green's functions, 



£&(*,*') = ^(*-0([c^W,4 ct (0]+)gc, 



(B2) 

and the occupation number of electrons in each quantum dot, 

(naa(t)) = (cL(*)<w(*)> GO - (B3) 

To do this, we attach fictitious field e~ lp °"'^ to for current, and add fictitious energy h a(T (t) to e aiJ for the average 
number. Then, from the generating functional W the above quantities can be calculated as, 



Iaait) = - iq 



Gaa(T,T') = 



5Ap aa (t) 



{n aa {t)) = -%h 



Sg ao -(T, t') 
SW 



5Ah m (t) 



Paa=h at7 =0 

) 

Pa<r=h aa =0 

(B4) 

Pa<r=h a cr=0 
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where p aa {t) — ±Ap„ rr ft)/2 and h arT (t) — ±Ah aa (t)/2 are assumed on the upper (+) and lower (— ) branches of 
the Keldysh contour. [22, [23| With the fictitious fields, the evaluation of the generating functional is straightforward 
because bare Green's function is simply changed as, 



0m 1 ( r > T ') -> (i^d T -e m [e aa +h atr (T)] -A)5(r-r')- 



(B5) 



In order to evaluate the functional derivatives, we expand the generating functional in series, W = X^^Lo 
where 



W (0) = ^ ?p Trln[ ff -V^], 



(B6) 



Here, <7 p (t, r') contains the fictitious fields, and Sp™^(r',r) represent all the proper and improper n-th order self- 
energies. By performing the functional derivatives we obtain, 



SW 



8Ah aa {t) 
SW 

SAp ai j (t) 



E 



de 



Sh aa (r) 
de aa J 5Ah aa (t) 

Sp aa (r) 



9m(T, t') + ,g m (r, (Ti,T2)g m (T2,T')\ drdr', 



5Ap acr (t) 



9acr(T, T')Y, a(T (T' , t) - S ao -(r, r')g aa (V, t) 



drdr'. 



(B7) 



By expressing t = J2„ ^ {n) = s + S .9 S + S.gS.gS + . . . = g~ l GJl = EGg" 1 with proper self-energy S, and performing 
the Keldysh rotation for the projection onto the real time, we finally obtain 

/oo 
dt' [G*(M')£^ (*',*) H-G^M')^ (*',*)] , 
-oo 

Gaa{T,T') = S Qcr (t, t' ) + / S Q(T (t, ti )G Qcr (n , t 2 ) (t 2 , t')<2ti dr 2 , 



(rw(i)) = iK/,< 



de m 
de ac , 



G<(t,t). 



(B8) 



For static cases, since Green's functions depend only on the difference between the time arguments, the expression 
for the current is further reduced in energy representation of Eq. (|32|) . 



/OO 
d£ [g> a (E)Z< a (E) - g< a (E)X> a (E)] , 



(B9) 



where in the second line we make use of Eq. (|30 



Appendix C: Current conservation 

The current conservation can be shown by concentrating on one of n-th order diagrams in the generating function, 
which consist of n conduction or 2n auxiliary particle Green's functions. Each diagram can be expressed either in 
terms of n conduction Green's functions, or in terms of 2n auxiliary particle Green's functions. Since they represent 
the same diagram, we can write 

Y,\ L a Ar,r')^(r',r)drdT' = - ]T gi j> 9m (r,r')t^ (t> ,r)drdr' . (CI) 
By summing all diagrams in the generating functional, we arrive at 

J2 fG aa (T,T')Z aa (T',T)dT> =-Y J <f>^G m {T,T , )Y, m {T , ,T)dT'. (C2) 



Additionally, by applying the Keldysh rotation onto real time, the above relation becomes 

m 

Using this relation, the sum of currents through both tunneling barriers can be written as 

£/ QCT (f) = -q^Y /^'[ G ™^OS<(^,0-G<(i,t')S>(t',t) 

a a m 

In static case, we obtain, by using the energy representation, 

^ Ia ° = ~2^ R £ Y I dE [ G ^( E )^{E) - G< (£)£>(£)] = 0, 

olo m 

where Eq. (|26|) is used. 



